Groupe: Équipe O
1 - Onorato, Claudia, Matricule: 1845448
2 - Harvey, William, Matricule: 1851388
Séance 1¶
Ce laboratoire a pour objectif de manipuler et analyser les signaux 1D dans le domaine temporel et fréquentiel.
Remise:
La date de remise est le vendredi 28 septembre à 23h59. Une pénalité de 3 points par jour sera appliquée lors d'un retard.
Documents à remettre :
Les exercices doivent être codés dans ce fichier TP.ipynb. Les réponses aux questions doivent être incluses dans le code sous forme de commentaires ou dans des cellules dédiées (Markdown ou text). Les exercices doivent être séparés par des cellules, suivant le template fourni. Vous devez bien identifier chaque exercice et sous-question, et bien commenter le code. Veuillez nommer vos variables de manière explicite et assurez-vous que toutes les figures soient lisibles.
Créer un fichier de rendu html (Fichier -> Télécharger au format... -> HTML) de votre code et de vos graphiques. Veuillez remettre tous vos fichiers (.ipynb, html et autres) dans un seul fichier zip et nommez ce fichier selon vos matricules (Mat1_Mat2.zip).
Une pénalité de 3 points sera appliquée si ces consignes ne sont pas respectées.
import numpy as np
import timeit
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (20,4)
Exercice 1 (2 points)¶
1-a. (0.25 point) Tracez le signal $s(t) = \dfrac{sin(\pi t)}{\pi t}$ sur l'intervalle $t \in [-4;4]$ avec une résolution de 1000 points (vous pouvez utiliser la fonction np.linspace pour générer le vecteur $t$). N'oubliez pas le titre et la légende (comme pour tous les autres graphiques que vous tracerez en laboratoire) !
t = np.linspace(-4, 4, num=1000)
s = np.sin(np.pi * t)/(np.pi * t)
plt.plot(t, s)
plt.grid()
plt.title("Amplitude d'une fonction sinus cardinal s(t) en fonction de t")
plt.xlabel('t')
plt.ylabel('s(t)')
plt.show()
1-b. (0.25 point) Comment se nomme ce signal ?
Le signal s(t) correspond à une fonction sinus cardinal.
- (0.5 point) Tracez dans une nouvelle figure la superposition des signaux: \begin{align*} s_1&=sin(t) &s_2&=\dfrac{\sin(3 t)}{3} & s_3&=\dfrac{\sin(5 t)}{5} \end{align*} sur l'intervalle $t \in [-2; 2]$ et avec une résolution de 500 points.
t = np.linspace(-2, 2, num=500)
s_1 = np.sin(t)
s_2 = np.sin(3 * t) / 3
s_3 = np.sin(5 * t) / 5
plt.plot(t, s_1, label='$s_1$')
plt.plot(t, s_2, label='$s_2$')
plt.plot(t, s_3, label='$s_3$')
plt.grid()
plt.title("Amplitude des signaux $s_x(t)$ en fonction de t")
plt.xlabel('t')
plt.ylabel('s(t)')
plt.legend()
plt.show()
- (0.5 point) Superposez à ces 3 signaux leur somme $s_1 + s_2 + s_3$. Pour rendre le diagramme plus lisible vous pouvez tracer les 3 signaux en pointillés en ajoutant l'argument '--' à la fonction
plt.plot(x, ..., '--'). Tracez la somme en gras (plt.plot(x, ..., linewidth=3)) sur le même graphique. Ajoutez une légende correspondant à chaque courbe.
plt.plot(t, s_1, '--', label='$s_1$')
plt.plot(t, s_2, '--', label='$s_2$')
plt.plot(t, s_3, '--', label='$s_3$')
plt.plot(t, (s_1 + s_2 + s_3), linewidth=3, label='somme')
plt.grid()
plt.title("Amplitude des signaux $s_x(t)$, et de leur somme, en fonction de t")
plt.xlabel('t')
plt.ylabel('s(t)')
plt.legend()
plt.show()
4-a. (0.25 point) Affichez dans une nouvelle figure mais sur le même intervalle $t\in[-2;2]$ et la même résolution de 500 points, le signal: $S_{50}(t) = \dfrac{1}{2} + \dfrac{2}{\pi} \sum_{i=0}^{50} \dfrac{\sin((2i+1) t)}{2i+1}$, c'est à dire la somme de $\frac{1}{2}$ et des $\frac{2}{\pi}\frac{\sin(k t)}{k}$ avec $k$ impair de $1$ à $101$
t = np.linspace(-2, 2, num=500)
def get_signal(n=50):
elements_before_sum = np.array([np.sin((2*i+1)*t)/(2*i+1) for i in range(n+1)])
summation = np.sum(elements_before_sum, axis=0)
return (1/2) + (2/np.pi) * summation
plt.plot(t, get_signal(), label='$s_50(t)$')
plt.grid()
plt.title("Amplitude du signal $s_50(t)$ en fonction de t")
plt.xlabel('t')
plt.ylabel('s(t)')
plt.legend()
plt.show()
4-b. (0.25 point) Enfin, réessayez pour $i$ allant de $0$ à $500$ ($k$ allant de $1$ à $1001$). Sur l'intervalle considéré, quel signal remarquable est approximé par cette somme de sinusoïdales ?
plt.plot(t, get_signal(n=500), label='s_500(t)')
plt.grid()
plt.title("Amplitude du signal $s_{500}(t)$ en fonction de t")
plt.xlabel('t')
plt.ylabel('s(t)')
plt.legend()
plt.show()
Réponse: Le signal remarquable approximé par cette somme de sinusoïdales est la fonction échelon. Celle-ci est décrite, par $s(t) = 1$ pour $t >= 0$, sinon $s(t) = 0$.
Exercice 2 (3.5 points)¶
Soit le signal analogique $Y(t)$ suivant :
\begin{align} Y(t) = 2\sin( 165\pi t ) + 13\cos( 6\pi t ) - 3\cos( 80\pi t ) \end{align}- (0.5 point) Déterminez théoriquement les fréquences présentes dans ce signal.
Réponse:
En utilisant: $$165 \pi t = 2 \pi ft$$ on peut isoler f: $$f = 82,5 Hz$$ On a donc la fréquence 82,5Hz qui compose notre signal.
Ensuite, une autre composante est additionnée à notre signal. $$6\pi t = 2 \pi ft$$ $$f = 3Hz$$
Finalement, on a aussi: $$80 \pi t = 2 \pi ft$$ $$f = 40Hz$$
En conclusion, notre signal est composé des fréquences 3, 40 et 82,5 Hz
- (1 point) Tracez le signal $Y(t)$ pour $0 \leq t \leq 1$ avec une fréquence d'échantillonnage $F_e$ de 20Hz, 75Hz, 100Hz, 160Hz, 180 Hz et 330Hz. Utilisez la commande plt.subplots afin de créer plusieurs graphique (partageant les mêmes abscisses) et de les afficher l'un en dessous de l'autre. Ajoutez un titre au graphique, une étiquette aux axes pour chaque graphique. Utilisez ylim pour contrôler l'échelle de l'ordonnée.
SAMPLING_FREQUENCIES = [20, 75, 100, 160, 180, 330]
T_MIN = 0 # in seconds
T_MAX = 1 # in seconds
plt.rcParams["figure.figsize"] = (20,15)
fig, axes = plt.subplots(nrows=len(SAMPLING_FREQUENCIES), sharex=True)
for sampling_frequency, axe in zip(SAMPLING_FREQUENCIES, axes):
t = np.linspace(T_MIN, T_MAX, num=sampling_frequency*(T_MAX-T_MIN))
y_t = 2 * np.sin(165*np.pi*t) + 13 * np.cos(6*np.pi*t) - 3 * np.cos(80*np.pi*t)
axe.plot(t, y_t, label=f'Y(t) échantilloné à {sampling_frequency}Hz')
axe.grid()
axe.set_ylabel('Y(t)')
axe.legend()
plt.xlabel('t (s)')
fig.suptitle("Signal $Y(t)$ selon différentes fréquences d'échantillonage")
fig.show()
- (1 point) Comment la fréquence d'échantillonnage affecte la forme du signal?
La fréquence d'échantillonage affecte la forme du signal tel que l'on ne peut pas, avec une basse fréquence d'échantillonage, voir les hautes fréquences dans le signal. En effet, si l'on regarde par exemple le signal $Y(t)$ échantilloné à 20 Hz, on ne peut distinguer qu'une seule fréquence (la plus basse).
- (1 point) Lesquelles, parmi ces fréquences d'échantillonnage, satisfont le théorème de Nyquist-Shannon? En pratique, quel compromis doit-on faire lors du choix d'une fréquence d'échantillonnage ?
Selon le théorème de Shannon, la fréquence d'échantillonage doit être au moins deux fois plus grande que la fréquence maximale qui compose le signal.
Dans le cas de $Y(t)$, la fréquence maximale $f_{m}$ est de 82,5 Hz. Ainsi, les fréquences d'échantillonage respectant le théorème de Nyquist-Shannon sont celles qui sont égales ou plus grandes à $2*82,5=165 Hz$. Les fréquences du numéro précédent qui respectent le théorème sont donc 180 et 330 Hz.
Donc, en pratique, le compromis lorsqu'on doit choisir une fréquence d'échantillonage est entre la précision du signal accumulé (qui augmente avec la fréquence d'échantillonage) et les contraintes en espaces mémoire et en temps de traitement. En effet, en augmentant la fréquence d'échantillonage, le signal recueilli sera plus volumineux, et si le signal doit être traité en temps réel, les calculs devront être plus optimisé afin de pouvoir suivre la cadence.
Exercice 3 (5.5 points)¶
Soient les trois signaux sinusoïdaux $Y_1(t)$, $Y_2(t)$ et $Y_3(t)$ suivants : \begin{align*} Y_1(t) &= 7 \sin( 2\pi\times10 t )\\ Y_2(t) &= 4 \sin( 2\pi\times25 t + \frac{\pi}{3})\\ Y_3(t) &= 3 \cos( 2\pi\times50 t ) \end{align*}
Ces signaux sont échantillonnés à la fréquence $F_e = 250$Hz et seront observés sur l'intervalle: \begin{align*} 0 \leq t \leq 1 \end{align*}
- (1 point) Tracez les trois signaux $Y_1(t)$, $Y_2(t)$ et $Y_3(t)$. Mettez un titre, une légende et les étiquettes des axes x et y.
plt.rcParams["figure.figsize"] = (20,5)
SAMPLING_FREQ = 250 # in Hz
T_MIN = 0 # in seconds
T_MAX = 1 # in seconds
t = np.linspace(T_MIN, T_MAX, num=SAMPLING_FREQ*(T_MAX-T_MIN))
y_1 = 7 * np.sin(2 * np.pi * 10 * t)
y_2 = 4 * np.sin(2 * np.pi * 25 * t + np.pi/3)
y_3 = 3 * np.cos(2 * np.pi * 50 * t)
plt.plot(t, y_1, label='$y_1$')
plt.plot(t, y_2, label='$y_2$')
plt.plot(t, y_3, label='$y_3$')
plt.grid(which='both')
plt.title("Amplitude des signaux $y_x(t)$ en fonction de t")
plt.xlabel('t')
plt.ylabel('$y_x$(t)')
plt.legend()
plt.show()
- (0.5 point) Déterminez graphiquement la période de chacun de ces signaux. Comparez chaque résultat avec sa valeur théorique.
En observant le graphique, on peut voir que
- la période de $y_1(t)$ est de $0.1 secondes$ (tel que $y_1(t+0,1)=y_1(t)$).
- la période de $y_2(t)$ est d'environ $0.04 secondes$ (tel que $y_2(t+0,04)=y_2(t)$).
- la période de $y_3(t)$ est d'environ $0.02 secondes$ (tel que $y_3(t+0,02)=y_3(t)$).
De manière analytique, on peut voir que:
- la période de $y_1(t)$ peut être déterminé en isolant $T$: $$2\pi*10t=2\pi(1/T)t$$ $$T=\frac{2\pi t}{2\pi10t}=1/10=0,1 Hz$$
- la période de $y_2(t)$ peut être déterminé en isolant $T$ (et en ignorant le décalage de la phase): $$T=\frac{2\pi t}{2\pi25t}=1/25=0,04 Hz$$
- la période de $y_2(t)$ peut être déterminé en isolant $T$: $$T=\frac{2\pi t}{2\pi50t}=1/50=0,02 Hz$$
3-a. (0.25 point) Tracez le signal composite $Z(t)=Y_1(t) + Y_2(t) + Y_3(t)$.
Z_t = y_1+y_2+y_3
plt.plot(t, Z_t, label='$Z(t)$')
plt.grid(which='both')
plt.title("Signal composite $Z(t)$ en fonction de t")
plt.xlabel('t')
plt.ylabel('Z(t)')
plt.legend()
plt.show()
3-b. (0.25 point) Graphiquement, quelle semble être la fréquence du signal $Z(t)$? Déterminez analytiquement cette fréquence sachant que la fréquence d'un signal composite est égal au plus grand diviseur commun des fréquences des signaux qui le composent.
Graphiquement, la période du signal correspond à 0,2, comme on peut voir que le signal se répète tel que $Z(t+0,2)=Z(t)$. Ainsi, la fréquence déduite est de $f=1/T=1/0,2=5 Hz$.
Sachant que les trois fréquences correspondantes aux trois signaux additionnés, qui constitue le signal $Z(t)$, sont $1/0,1=10 Hz$, $1/0,04=25 Hz$ et $1/0,02=50 Hz$, on peut déduire la fréquence du signal composite. En effet, le plus grand commun diviseur est bien $5 Hz$, ce qui correspond à la fréquence qui a été déduite graphiquement.
- (1 point) Calculez la transformée de Fourier rapide (FFT) des signaux $Y_1(t)$, $Y_2(t)$ et $Y_3(t)$ à l'aide de la fonction np.fft.fft et affichez le spectre de fréquence de chacun. Que remarquez vous?
freqs = np.fft.fftfreq(t.shape[-1], d=(1/SAMPLING_FREQ))
xticks = np.arange(np.min(freqs), np.max(freqs), 5)
for idx, y_x in enumerate([y_1,y_2,y_3]):
plt.plot(freqs, np.fft.fft(y_x))
plt.xticks(xticks)
plt.grid(which='both')
plt.title(f"Spectre du signal $Y_{idx+1}(t)$")
plt.xlabel('f (Hz)')
plt.ylabel(f'TFD[$Y_{idx+1}(t)$]')
plt.show()
- (1 point) Calculez la TFD du signal composite et affichez son spectre de fréquence. Que remarquez-vous?
plt.plot(freqs, np.fft.fft(Z_t))
plt.xticks(xticks)
plt.grid(which='both')
plt.title(f"Spectre du signal composite $Z(t)$")
plt.xlabel('f (Hz)')
plt.ylabel('TFD[Z(t)]')
plt.show()
# On remarque que le spectre du signal composite Z(t) correspond à la somme des spectres des signaux qui le compose.
(1 point) On rappelle que la transformée de Fourier discrète (TFD) du signal s(n) à N échantillons est définie par: \begin{align} S[k] = \sum_{n=0}^{N-1} s(n)e^{-2i\pi n k/N} \end{align} pour $0 \leq k < N$. Jusqu'à présent, vous avez utilisé la fonction fft (Fast Fourier Transform) offerte par numpy, qui n'est qu'une optimisation computationnelle du calcul de la TFD. Mais est-elle vraiment plus rapide que la TFD brute? Complétez le code suivant calculant la TFD d'un signal s:
def tfd(s): N = len(s) S = np.zeros(N, dtype=complex) n = np.arange(0, N) exp_basis = np.exp(-2j*np.pi*n/N) for k in range(N): S[k] = np.sum(...) return S
Vérifiez que vous obtenez le même spectre de $Z(t)$ avec votre fonction et la fft (tracez-le).
def tfd(s):
N = len(s)
S = np.zeros(N, dtype=complex)
n = np.arange(0, N)
exp_basis = np.exp(-2j*np.pi*n/N)
for k in range(N):
S[k] = np.sum(s*np.power(exp_basis, k))
return S
plt.plot(freqs, tfd(Z_t))
plt.xticks(xticks)
plt.grid(which='both')
plt.title(f"Spectre du signal $Z(t)$ calculé manuellement")
plt.xlabel('f (Hz)')
plt.ylabel('TFD[Z(t)]')
plt.show()
- (0.5 point) A l'aide du module timeit, comparez le temps d'executions de la fonction np.fft.fft et de votre fonction tdf, sur 100 itérations. Aidez-vous pour cela du code ci-dessous. Qu'en concluez-vous?
temps = timeit.timeit(lambda:ma_fonction(args), n=nb_iterations) # A ajuster
## Code ##
NB_ITERATION = 100
temps_tfd = timeit.timeit(lambda:tfd(Z_t), number=NB_ITERATION)
temps_fft = timeit.timeit(lambda:np.fft.fft(Z_t), number=NB_ITERATION)
print(f"Temps d'exécution de notre fonction TDF sur {NB_ITERATION} itérations: {temps_tfd}")
print(f"Temps d'exécution de la fonction np FFT sur {NB_ITERATION} itérations: {temps_fft}")
Réponse: Nous pouvons en conclure qu'il est beaucoup plus rapide d'utiliser la fast fourier transform comparativement à la transformée de Fourier discrète.
Un mystérieux enregistrement a été retrouvé lors d'une maintenance des serveurs Moodle. Cependant il semble bien peu audible; comme si un chargé de laboratoire facétieux s'était amusé à y ajouter des signaux parasites...
Qu'à cela ne tienne, vous vous proposez de le décrypter en utilisant votre expertise en traitement du signal.
import IPython.display as ipd
import scipy.io.wavfile as wavfile
from scipy.signal import firwin, lfilter, freqz
- (1 point) Chargez le fichier audio.wav à l'aide de la fonction scipy.io.wavfile.read. Celle-ci renvoie deux variables: la fréquence d'échantillonnage et le signal audio.
Vous pouvez écouter le signal sonore original ou filtré à l'aide de la fonction IPython.display.Audio.
ipd.Audio(signal, rate=sample_rate)
En écoutant le signal, repérez la ou les perturbations de la mélodie principale. Pour chaque perturbation, précisez si le signal est haute ou basse fréquence.
FILE_NAME = 'corrupted_audio.wav'
sample_rate, signal = wavfile.read(FILE_NAME)
print(f"Sample rate is {sample_rate}Hz and audio lasts {signal.shape[-1]/sample_rate} seconds")
ipd.Audio(signal, rate=sample_rate)
Réponse:
Il semble y avoir deux signaux parasites de fréquence différente. Ceux-ci correspondent à deux son strident, un plutôt aïgu (haute fréquence) et l'autre plus grave (basse fréquence).
2-a. (0.5 point) Calculez la TFD du signal audio et affichez-la. Ajouter un titre et les axes.
freqs = np.fft.fftfreq(signal.shape[-1], d=(1/sample_rate))
xticks = np.arange(np.min(freqs), np.max(freqs), 1500)
plt.plot(freqs, np.fft.fft(signal))
plt.xticks(xticks)
# plt.xlim(-2400,2400) # helps to see what are the pertubations' frequencies
plt.grid(which='both')
plt.title(f"Spectre du signal audio")
plt.xlabel('f (Hz)')
plt.ylabel('FFT')
plt.show()
2-b. (0.5 point) Quelle est la note (do, ré, mi, fa, sol, la ou si) correspondant aux perturbations? Aidez-vous de l'article Wikipédia et du spectre de fréquence affiché.
Les notes correspondantes aux pertubations sont le Ré de la 6e octave pour la perturbation à la fréquence 2 349,32 Hz et le Do de la 3e octave pour la perturbation ayant la fréquence 261,63 Hz.
RÉ_OCTAVE_6 = 2349.32
DO_OCTAVE_3 = 261.63
plt.plot(freqs, np.fft.fft(signal))
plt.axvline(RÉ_OCTAVE_6, c='r', ls='--')
plt.axvline(-RÉ_OCTAVE_6, c='r', ls='--')
plt.axvline(DO_OCTAVE_3, c='r', ls='--')
plt.axvline(-DO_OCTAVE_3, c='r', ls='--')
plt.xticks(xticks)
# plt.xlim(0,2400) # helps to see what are the pertubations' frequencies
plt.grid(which='both')
plt.title(f"Spectre du signal audio")
plt.xlabel('f (Hz)')
plt.ylabel('FFT[Z(t)]')
plt.show()
- (0.5 point) Pour filtrer la perturbation la plus aigue, vous allez concevoir un filtre passe-bas. Comme vu en cours, un filtre est défini par plusieurs paramètres: sa fréquence de coupure, son ordre et le type de fenêtre employé. Donnez la définition de la fréquence de coupure.
Fréquence de coupure pour un filtre idéal: Dans le cas d'un filtre passe-bas, la fréquence de coupure est la fréquence à partir de laquelle le gain du filtre passe de 1 à 0. Pour un filtre passe-haut, c'est l'inverse. La fréquence de coupure est la fréquence à partir de laquelle le gain du filtre passe de 0 à 1.
Fréquence de coupure pour un filtre réel: En réalité, un filtre ne peut pas directement passer de 0 à 1 instantannément. Ainsi, dans le cas des filtres réels, la fréquence de coupure est la fréquence à partir de laquelle le gain est de 0.5 (-3 dB).
- (1 point) En utilisant les fonctions scipy.signal.firwin et scipy.signal.lfilter, concevez un filtre passe bas d'ordre 128 pour filtrer la perturbation à la plus haute fréquence. En tenant compte de la définition donnée ci-dessus, comment choisir une bonne fréquence de coupure? Filtrez la musique avec le filtre que vous venez de concevoir et décrivez ce que vous entendez. N'hésitez pas à jouer sur la fréquence de coupure pour observer son effet. Que se passe-t-il si elle est trop haute? Et si elle est trop basse? Quel type de filtre serait plus approprié à la pertubation considérée?
FILTER_ORDER = 128
CUTOFF_FREQ = RÉ_OCTAVE_6 - 1000
fir_filter = firwin(
numtaps=FILTER_ORDER+1, # Documentation: numtaps equals to the filter order + 1
cutoff=CUTOFF_FREQ,
fs=sample_rate,
pass_zero='lowpass'
)
filtered_signal_window_method = lfilter(fir_filter, [1.0], signal)
ipd.Audio(filtered_signal_window_method, rate=sample_rate)
# time-domain filter kernel
plt.plot(fir_filter)
plt.xlabel('Time points')
plt.title('Filter kernel (firwin)')
plt.show()
# frequency-domain filter kernel
kernel_freqs = np.fft.fftfreq(fir_filter.shape[-1], d=(1/sample_rate))
plt.plot(kernel_freqs, np.abs(np.fft.fft(fir_filter))**2, 'bo-', label='Actual')
plt.plot([np.min(kernel_freqs),-CUTOFF_FREQ,-CUTOFF_FREQ,CUTOFF_FREQ,CUTOFF_FREQ,np.max(kernel_freqs)],[0,0,1,1,0,0],'ro-',label='Ideal')
plt.grid(which='both')
plt.title("Spectre du filtre ayant $F_{coupure}$= " + f"{CUTOFF_FREQ:.2f}")
plt.xlabel('f (Hz)')
plt.ylabel('TFD')
plt.legend()
plt.show()
# frequency-domain filtered signal
plt.plot(freqs[0:freqs.shape[0]//2], np.fft.fft(filtered_signal_window_method)[0:freqs.shape[0]//2])
plt.grid(which='both')
plt.xlim(0,RÉ_OCTAVE_6+20) # helps to see what are the pertubations' frequencies
# plt.ylim(-3000,3000) # helps to see what are the pertubations' frequencies
plt.title("Spectre du signal audio filtré ayant $F_{coupure}$= " + f"{CUTOFF_FREQ:.2f}")
plt.xlabel('f (Hz)')
plt.ylabel('TFD')
plt.show()
Réponse:
En considérant que la fréquence de coupure est la fréquence pour laquelle le gain est de 0,5, on peut choisir une fréquence de coupure, dans le cas d'un filtre passe-bas, un peu plus basse que la fréquence du signal parasite. Ainsi, on veut qu'à partir de la fréquence parasite, le gain soit de 0, et qu'avant celle-ci, elle ne soit pas à 0.
Si on fixe la fréquence de coupure à une valeur trop basse (soit beaucoup plus basse que $2349.32 Hz$), on n'entend plus la fréquence aigüe, mais on entend également moins la musique que l'on voudrait entendre.
Si on fixe la fréquence de coupure à une valeur trop basse (soit environ $2349.32 Hz$ ou plus), on ne filtre pas ou on filtre un peu la fréquence aigüe. On entend alors encore celle-ci, ou elle est un peu moins bruyante.
Un filtre qui serait plus approprié pour enlever cette fréquence parasite du signal serait un notch filter, qui est un coupe-bande sur un intervalle très étroit. En effet, en zoomant sur le spectre du signal, on voit qu'il se trouve des amplitudes non nulles aux fréquences au-delà de la fréquence du signal parasite. En faisant un passe-bas, on enlève ces fréquences du signal. On perd alors des fréquences du signal d'intérêt. En faisant un notch filter, on peut enlever une fréquence spécifique, comme son intervalle est très court.
- (1 point) Créez trois filtres passe-haut d'ordre 128 et de fréquence de coupure de 750 Hz en utilisant les fenêtres suivantes : Chebyshev, Hamming et Blackman. Pour la fenêtre Chebyshev, utilisez une atténuation de 30 dB.
## Code ##
FILTER_ORDER = 128
CUTOFF_FREQ = 750 # Hz
CHEBY_ATTENUATION = 30 # dB
hamming_fir_filter = firwin(
numtaps=FILTER_ORDER+1,
cutoff=CUTOFF_FREQ,
window='hamming',
fs=sample_rate,
pass_zero='highpass'
)
blackman_fir_filter = firwin(
numtaps=FILTER_ORDER+1,
cutoff=CUTOFF_FREQ,
window='blackman',
fs=sample_rate,
pass_zero='highpass'
)
cheby_fir_filter = firwin(
numtaps=FILTER_ORDER+1,
cutoff=CUTOFF_FREQ,
window=('chebwin', CHEBY_ATTENUATION),
fs=sample_rate,
pass_zero='highpass'
)
- (1.5 points) Affichez, pour chaque filtre, sa réponse fréquentielle à l'aide de la fonction scipy.signal.freqz. Sur chaque graphique, tracez la droite verticale passant par la fréquence de coupure. Vous pouvez pour cela utiliser la fonction:À quelle atténuation coupe-t-elle la courbe de la réponse fréquentielle? Était-ce prévisible?
plt.axvline(freq_coupure, color='r', linewidth=2)
for win_name, win_filter in [
('hamming', hamming_fir_filter),
('blackman', blackman_fir_filter),
('chebyshev', cheby_fir_filter)
]:
w,h = freqz(win_filter, fs=sample_rate)
plt.plot(w, np.abs(h) , 'b')
plt.axvline(CUTOFF_FREQ, color='r', linewidth=2)
plt.axhline(0.5, color='g',ls='--') # added gain at frequency cutoff [as answered at question 3]
plt.grid(which='both')
plt.title(f"FIR filter with {win_name} window frequency response")
plt.xlabel('Frequence (Hz)')
plt.ylabel('Gain')
plt.show()
Réponse:
Comme on peut voir ci-haut à l'aide de la ligne verte, qui se trouve à $0,5$, la fréquence de coupure, pour les trois filtres, correspond à une atténuation de $0,5$. On pouvait s'y attendre, comme on a écrit plus haut à la question 3, car la définition de la fréquence de coupure est au point où le signal a un gain de 50%, ce qui correspond à -3 dB.
- (1.5 points) Filtrez le signal audio qui a déjà été filtré avec le filtre passe-bas à l'aide des trois filtres créés précédemment. Lorsque vous écoutez les trois signaux, que remarquez-vous? Pour le filtre passe-haut et en fonction de la fréquence de coupure utilisée, quel va être le compromis sur la qualité du signal restaurée ?
for win_name, win_filter in [
('hamming', hamming_fir_filter),
('blackman', blackman_fir_filter),
('chebyshev', cheby_fir_filter)
]:
currently_filtered_signal = lfilter(win_filter, [1.0], filtered_signal_window_method)
print(f'Audio sample filtered with window {win_name}')
ipd.display(ipd.Audio(currently_filtered_signal, rate=sample_rate))
Réponse:
Nous pouvons remarquer que les trois signaux filtrés possèdent encore les signaux parasites. En effet, le signal parasite à la fréquence basse est encore très audible. Par contre, comparativement au signal audio résultant de l'application du premier filtre passe-bas, on entend davantage la musique en arrière-plan. On peut se douter que l'on entend encore la basse fréquence par le fait que la fréquence de coupure est trop basse; en effet, si l'on augmente la fréquence de coupure à 1500 Hz, le gain à 261.63 Hz est environ de 0, et la fréquence basse n'est presque plus audible. On voit également que le signal résultant est différent dépendamment de la fenêtre choisie.
Le compromis est encore une fois au niveau des fréquences du signal orignal qui seront perdues lors de l'application du filtre. En effet, on réussit à beaucoup plus réduire la fréquence basse parasite à 261.63 Hz en choisissant une fréquence de coupure à 1500 Hz plutôt qu'à 750 Hz, mais on perd une grande partie des fréquences du signal d'intérêt.
- (1 point) Calculez les TFD des signaux filtrés avec les trois filtres et affichez les spectres. Est-ce que les spectres, au niveau des basses fréquences, correspondent à ce que vous avez entendu?
for win_name, win_filter in [
('hamming', hamming_fir_filter),
('blackman', blackman_fir_filter),
('chebyshev', cheby_fir_filter)
]:
currently_filtered_signal = lfilter(win_filter, [1.0], filtered_signal_window_method)
plt.plot(freqs[0:freqs.shape[0]//2], (np.fft.fft(currently_filtered_signal))[0:freqs.shape[0]//2])
plt.grid(which='both')
plt.xlim(0, np.max(freqs))
plt.title(f"Spectre du signal audio filtré avec {win_name}")
plt.xlabel('f (Hz)')
plt.ylabel('TFD')
plt.show()
# pour mieux observer les basses fréquences, on réaffiche les spectres en limitant l'axe des fréquences
for win_name, win_filter in [
('hamming', hamming_fir_filter),
('blackman', blackman_fir_filter),
('chebyshev', cheby_fir_filter)
]:
currently_filtered_signal = lfilter(win_filter, [1.0], filtered_signal_window_method)
plt.plot(freqs[0:freqs.shape[0]//2], (np.fft.fft(currently_filtered_signal))[0:freqs.shape[0]//2])
plt.grid(which='both')
plt.xlim(0, RÉ_OCTAVE_6+20) # pour mieux observer les basses fréquences
plt.title(f"Spectre du signal audio filtré avec {win_name}")
plt.xlabel('f (Hz)')
plt.ylabel('TFD')
plt.show()
for win_name, win_filter in [
('hamming', hamming_fir_filter),
('blackman', blackman_fir_filter),
('chebyshev', cheby_fir_filter)
]:
currently_filtered_signal = lfilter(win_filter, [1.0], filtered_signal_window_method)
plt.plot(
freqs[0:freqs.shape[0]//2],
10*np.log10((np.fft.fft(currently_filtered_signal))[0:freqs.shape[0]//2]),
label=f"Spectre du signal audio filtré avec {win_name}"
)
plt.title('Spectre en échelle logarithme des signaux filtrés')
plt.legend()
plt.grid(which='both')
plt.xlim(0, np.max(freqs))
# plt.xlim(0, RÉ_OCTAVE_6+20) # pour mieux observer les basses fréquences
plt.xlabel('f (Hz)')
plt.ylabel('TFD [dB]')
plt.show()
Réponse:
Afin de mieux répondre à la question, nous avons limité, dans des seconds graphiques, l'axe des fréquences qu'aux basses fréquences. Nous avons également superposé les spectres en échelle logarithmique pour mieux les comparer.
Oui, les spectres correspondent bien à quoi nous nous attendions. Nous étions d'abord surpris de voir des pics aux fréquences parasites, mais en comparant l'échelle à celle du spectre du signal filtré avec le passe-bas uniquement, on voit que les fréquences parasites ont été beaucoup atténuées. De plus, on peut également confirmer notre intuition que l'on entend plus clairement le signal d'intérêt, comme leurs amplitudes ne sont plus si petites comparativement aux fréquences parasites. En effet, si l'on regarde le signal filtré avec la fenêtre Chebyshev, on voit que les amplitudes du signal sont relativement grandes, ce qui transparaît dans l'extrait audio, comme c'est dans celui-là que l'on entend mieux la musique de fond.